#include "forcasIndependentes.h"

ForcasIndependentes::ForcasIndependentes(MODEL * md)
{
    pModel=md;
    MaterialOceanico = (md->SeaMat != NULL);
    if (md->SeaMat != NULL) DensidadeMatOc=pModel->SeaMat->Density;
    gravidade[0]=pModel->Gravity[0];
    gravidade[1]=pModel->Gravity[1];
}
ForcasIndependentesParalelo::ForcasIndependentesParalelo(MODEL * md)
{
    pModel=md;
    MaterialOceanico = (md->SeaMat != NULL);
    if (md->SeaMat != NULL) DensidadeMatOc=pModel->SeaMat->Density;
    gravidade[0]=pModel->Gravity[0];
    gravidade[1]=pModel->Gravity[1];
}

void ForcasIndependentes::Empuxo(SPHERE*s)
{
    if (MaterialOceanico)
    {
        REAL invDens;
        REAL raioQuad;
        REAL normaV;
        REAL mult;
        invDens = -1.0 * DensidadeMatOc;
        raioQuad = s->Rad * s->Rad;
        normaV = NORM(s->V);
        mult = invDens*raioQuad;
        mult *= normaV;
        s->Force( mult, s->V );// Equação do arrasto (Lord-Rayleigh)
    }
}
void ForcasIndependentes::Gravidade(SPHERE*s)
{
    s->Force( s->Mass, gravidade );
}

void ForcasIndependentes::Explosao(SPHERE*s)
{
	// Constantes
	REAL f = 1.7;
	REAL beta = 0.47;
	REAL rho = 9935;
	REAL n = 2.75;
	REAL c = 500;
	REAL W = 1.0;  // carga equivalente
	REAL R;   // distancia da particula a origem (onde esta a carga)
	REAL t;  // tempo
	REAL forcaP; // forca devido a pressao

	REAL Pos[2];
    Pos[0]=0.01-s->R[0];
    Pos[1]=0.01-s->R[1];
	R=NORM(Pos);
	Pos[0]/=R;
	Pos[1]/=R;
    forcaP = ( ( -f * beta * rho * pow(R/pow(W,0.3333),-n) * n*exp(-t/R*c)) / R + f * beta * rho * pow(R/pow(W,0.3333),-n) * t * c * exp(-t/R*c) /(R*R) ) * CONST_PI * s->Rad * s->Rad;
    s-> Force( forcaP, Pos );
}

void ForcasIndependentes::calcular()
{
    #pragma omp single
    {
        for (SPHERE* s = pModel->SPHList.Begin; s != NULL;s = s->Next)
        {
            Gravidade(s);
            //Explosao(s);
            Empuxo(s);
        }
    }
    #pragma omp barrier
}

void ForcasIndependentesParalelo::calcular()
{
    #pragma omp for schedule (dynamic,1)
    for (int i = 0; i < pModel->numProceEsf; i++)
    {
        SPHERE* testeFim;
        if (i!=pModel->numProceEsf-1) testeFim = pModel->threads[i+1];
        else testeFim=NULL;
        for (SPHERE* s = pModel->threads[i]; s != testeFim; s = s->Next)
        {
            Gravidade(s);
            //Explosao(s);
            Empuxo(s);
        }
    }
    //barreira implícita
}
